{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Unbinned Likelihood Tutorial\n",
    "\n",
    "The detection, flux determination, and spectral modeling of Fermi LAT sources is accomplished by a maximum likelihood optimization technique as described in the [Cicerone](https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/Cicerone_Likelihood/) (see also e.g. [Abdo, A. A. et al. 2009, ApJS, 183, 46](http://arxiv.org/abs/0902.1340)).\n",
    "\n",
    "To illustrate how to use the Likelihood software, this narrative gives a step-by-step description for performing an unbinned likelihood analysis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Unbinned vs Binned Likelihood\n",
    "\n",
    "**Unbinned** likelihood analysis is the preferred method for time series analysis of the LAT data, where the number of events in each time bin is expected to be small.\n",
    "\n",
    "However, for large time bins, analysis which includes bright background sources (such as the Galactic plane), and long time-baseline spectral and spatial analyses, a **binned** analysis is recommended.\n",
    "\n",
    "To perform a binned likelihood analysis, see the [Binned Likelihood](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/binned_likelihood_tutorial.html) tutorial. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Additional references**:\n",
    "* [SciTools References](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/references.html)\n",
    "* Descriptions of available [Spectral and Spatial Models](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/source_models.html)\n",
    "* Examples of [XML Model Definitions for Likelihood](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/xml_model_defs.html#xmlModelDefinitions):\n",
    "    * [Power Law](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/xml_model_defs.html#powerlaw)\n",
    "    * [Broken Power Law](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/xml_model_defs.html#brokenPowerLaw)\n",
    "    * [Broken Power Law 2](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/xml_model_defs.html#powerLaw2)\n",
    "    * [Log Parabola](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/xml_model_defs.html#logParabola)\n",
    "    * [Exponential Cutoff](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/xml_model_defs.html#expCutoff)\n",
    "    * [BPL Exponential Cutoff](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/xml_model_defs.html#bplExpCutoff)\n",
    "    * [Gaussian](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/xml_model_defs.html#gaussian)\n",
    "    * [Constant Value](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/xml_model_defs.html#constantValue)\n",
    "    * [File Function](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/xml_model_defs.html#fileFunction)\n",
    "    * [Band Function](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/xml_model_defs.html#bandFunction)\n",
    "    * [PL Super Exponential Cutoff](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/xml_model_defs.html#plSuperExpCutoff)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Prerequisites\n",
    "\n",
    "You will need an **event** data file, a **spacecraft** data file (also referred to as the \"pointing and livetime history\" file), and the current **background models** (available for [download](https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html)).\n",
    "\n",
    "You may choose to select your own data files, or to use the files provided within this tutorial.\n",
    "\n",
    "Custom data sets may be retrieved from the [LAT Data Server](http://fermi.gsfc.nasa.gov/cgi-bin/ssc/LAT/LATDataQuery.cgi)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Steps\n",
    "\n",
    "1. **Make Subselections from the Event Data**\n",
    "\n",
    "    Since there is computational overhead for each event associated with each diffuse component, it is useful to filter out any events that are not within the extraction region used for the analysis.\n",
    "\n",
    "\n",
    "2. **Make Counts Maps from the Event Files**\n",
    "    \n",
    "    By making simple FITS images, we can inspect our data and pick out obvious sources.\n",
    "\n",
    "\n",
    "3. **Download the latest diffuse models**\n",
    "\n",
    "    The recommended models for a normal point source analysis are `gll_iem_v07.fits` (a very large file) and `iso_P8R3_SOURCE_V2_v1.txt`. All of the background models along with a description of the models are available [here](https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html).\n",
    "\n",
    "\n",
    "4. **Create a Source Model XML File**\n",
    "\n",
    "    The source model XML file contains the various sources and their model parameters to be fit using the **gtlike** tool.\n",
    "    \n",
    "\n",
    "5. **Compute the diffuse responses**\n",
    "\n",
    "    Each event must have a separate response precomputed for each diffuse component in the source model.\n",
    "\n",
    "\n",
    "6. **Make an Exposure Map**\n",
    "\n",
    "    This is needed in order to analyze diffuse sources and derive absolute fluxes.\n",
    "    \n",
    "\n",
    "7. **Perform the Likelihood Fit**\n",
    "\n",
    "    Fitting the data to the model provides flux, errors, spectral indices, and other information.\n",
    "\n",
    "\n",
    "8. **Make Test-Statistic Maps**\n",
    "\n",
    "    These are used for point source localization and for finding weaker sources after the stronger sources have been modeled."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Make Subselections from the event data\n",
    "\n",
    " For this case the original extraction of data from the first six months of the mission was done as described in the Extract LAT Data tutorial.\n",
    "\n",
    "Selection of data:\n",
    "\n",
    "    Search Center (RA, DEC) =(193.98, -5.82)\n",
    "    Radius = 20 degrees\n",
    "    Start Time (MET) = 239557417 seconds (2008-08-04 T15:43:37)\n",
    "    Stop Time (MET) = 255398400 seconds (2009-02-04 T00:00:00)\n",
    "    Minimum Energy = 100 MeV\n",
    "    Maximum Energy = 500000 MeV\n",
    "\n",
    " We provide the user with the original photon and spacecraft data files extracted from the LAT Data Server:\n",
    "\n",
    "    L181204121625F588C43960_PH00.fits\n",
    "    L181204121625F588C43960_PH01.fits\n",
    "    L181204121625F588C43960_SC00.fits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/dataPreparation/L181204121625F588C43960_PH00.fits\n",
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/dataPreparation/L181204121625F588C43960_PH01.fits\n",
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/dataPreparation/L181204121625F588C43960_SC00.fits"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you need to combine multiple events files for your analysis, you must first generate a text file listing the events files to be included."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!mkdir data\n",
    "!mv *_PH* *_SC00* ./data\n",
    "!mv ./data/*_SC00* ./data/spacecraft.fits\n",
    "!ls ./data/*_PH* > ./data/events.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cat ./data/events.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This text file (`events.txt`) will be used in place of the input fits filename when running **gtselect**. The syntax requires that you use an @ before the filename to indicate that this is a text file input rather than a fits file.\n",
    "\n",
    "For simplicity, we rename the spacecraft file (the file ending in `SC00`) to `spacecraft.fits`.\n",
    "\n",
    "When analyzing a point source, it is recommended that you include events with high probability of being photons. To do this, you should use gtselect to cut on the event class, keeping only the SOURCE class events (event class 128) (or as recommended in the Cicerone).\n",
    "\n",
    "In addition, the full set of events have been divided into event types, which allow us to select events based on the conversion type, the quality of the track reconstruction, or the quality of the energy measurement. For standard selections, including all the event types for a given class, use `evtype=3` (\"INDEF\" is the default in **gtselect**)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```bash\n",
    "gtselect evclass=128 evtype=3\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " Be aware that `evclass` and `evtype` are hidden parameters. So to use them, you must type them on the command line.\n",
    "\n",
    "We perform a selection to filter the data we want to analyze. For this example, we consider all SOURCE class events (i.e. `evclass=128` and `evtype=INDEF`) within a 20 degree region of interest (ROI) centered on the blazar 3C 279. We apply the gtselect tool to the data file as follows: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtselect evclass=128 evtype=3\n",
    "    @./data/events.txt\n",
    "    ./data/3C279_region_filtered.fits\n",
    "    193.98\n",
    "    -5.82\n",
    "    20\n",
    "    INDEF\n",
    "    INDEF\n",
    "    100\n",
    "    500000\n",
    "    90"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " In the last step we also selected the energy range and the maximum zenith angle value (90 degrees). The Earth's limb is a strong source of background gamma rays. We filter them out with a zenith-angle cut. The value of 90 degrees is the one recommended by the LAT instrument team for analysis above 100 MeV. The filtered data are provided here.\n",
    "\n",
    "After the data selection is made, we need to select the good time intervals in which the satellite was working in standard data taking mode and the data quality was good. For this task we use gtmktime to select GTIs by filtering on information provided in the spacecraft file. The current gtmktime filter expression recommended by the LAT team in the Cicerone is: "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "(DATA_QUAL>0)&&(LAT_CONFIG==1)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " This excludes time periods when some spacecraft event has affected the quality of the data, ensures the LAT instrument was in normal science data-taking mode.\n",
    "\n",
    "Here is an example of running [gtmktime](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtmktime.txt) for our analysis of the region surrounding 3C 279."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtmktime\n",
    "    ./data/spacecraft.fits\n",
    "    (DATA_QUAL>0)&&(LAT_CONFIG==1)\n",
    "    no\n",
    "    ./data/3C279_region_filtered.fits\n",
    "    ./data/3C279_region_filtered_gti.fits"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " The data with all the cuts described above is provided in this link. A more detailed discussion of data selection can be found in the Data Preparation analysis thread.\n",
    "\n",
    "To view the DSS keywords in a given extension of a data file, use the gtvcut tool and review the data cuts on the EVENTS extension. This provides a listing of the keywords reflecting each cut applied to the data file and their values, including the entire list of GTIs. Here we use the option `suppress_gtis=yes` to keep the output easily readable: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtvcut suppress_gtis=yes\n",
    "    ./data/3C279_region_filtered_gti.fits\n",
    "    EVENTS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here you can see the event class selection, the location and radius of the data selection, the event type selection, as well as the energy range in MeV and the zenith angle cut. The time cuts to be used in the exposure calculation are defined by the GTI table.\n",
    "\n",
    "Various Fermitools will be unable to run if you have multiple copies of a particular DSS keyword. This can happen if the position used in extracting the data from the data server is different than the position used with [gtselect](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtselect.txt). It is wise to review the keywords for duplicates before proceeding. If you do have keyword duplication, it is advisable to regenerate the data file with consistent cuts.\n",
    "\n",
    "It is a good practice to check the maximum energy of the photons contained in our data sample. In this particular case, since we applied a maximum energy cut of 500 GeV, our statistics should be fine. However, if we try to go to higher energies, we might run out of photons and the likelihood fit will fail. To check the maximum energy, we can open the event file using e.g. *fv*:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!fv ./data/3C279_region_filtered_gti.fits"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Click on the `Hist` button for the EVENTS extension. Then select for the X axis `ENERGY`.\n",
    "\n",
    "The minimum and maximum energies are given by `Data Min` and `Data Max`. If Data Max is of the order of the maximum energy we used in our data selection we can proceed.\n",
    "\n",
    "However, if Data Max is much smaller than the energy used in **gtselect**, then we need to come back to **gtselect** and make a tighter energy cut to make sure we will have enough statistics. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/images/Likelihood/energy_max.png'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Make a counts map from the event data\n",
    "\n",
    "Next, we create a counts map of the ROI, summed over photon energies, in order to identify candidate sources and to ensure that the field looks sensible as a simple sanity check. For creating the counts map, we will use the gtbin tool with the option `CMAP` as shown below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtbin\n",
    "    CMAP\n",
    "    ./data/3C279_region_filtered_gti.fits\n",
    "    ./data/3C279_region_filtered_gti_cmap.fits\n",
    "    NONE\n",
    "    160\n",
    "    160\n",
    "    0.25\n",
    "    CEL\n",
    "    193.98\n",
    "    -5.82\n",
    "    0.0\n",
    "    AIT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">**NOTE**: Here [gtbin](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtbin.txt) shows a warning due to the fact we are not using the spacecraft file. This is not a concern.\n",
    "\n",
    "Now we can use the *ds9* command to launch a visualization tool:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!ds9 ./data/3C279_region_filtered_gti_cmap.fits"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*ds9* produces this display:\n",
    "\n",
    "<img src='https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/images/Likelihood/3C279_counts_map.png'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see one strong source and several weaker sources in this map. Mousing over the positions of these sources shows that two of them are likely 3C 279 and 3C 273.\n",
    "\n",
    "It is important to inspect your data prior to proceeding to verify that the contents are as you expect. A malformed data query or improper data selection can generate a non-circular region, or a file with zero events.\n",
    "\n",
    "By inspecting your data prior to analysis, you have an opportunity to detect such issues early in the analysis. (A more detailed discussion of data exploration can be found in the [Explore LAT Data](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/explore_latdata.html) analysis thread.) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Make an exposure map\n",
    "\n",
    " We are now ready to create an exposure map. The type of exposure map used by **Likelihood** differs significantly from the usual notion of exposure maps, which are essentially integrals of effective area over time. The exposure calculation that **Likelihood** uses consists of an integral of the total response over the entire ROI. See more information in the Cicerone.\n",
    "\n",
    "Since the exposure calculation involves an integral over the ROI, separate exposure maps must be made for every distinct set of DSS cuts. This is important if, for example, one wants to subdivide an observation to look for secular flux variations from a particular source or sources. To view the DSS keywords in a given extension of a data file, use the [gtvcut](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtvcut.txt) tool and review the data cuts for the EVENTS extension. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### a. Generate a livetime cube\n",
    "\n",
    " There are two tools needed for generating exposure maps. The first is gtltcube. This tool creates a livetimecube, which is a HealPix table, covering the full sky, of the integrated livetime as a function of inclination with respect to the LAT z-axis.\n",
    " \n",
    " For more information about the livetime cubes, see the documentation in the [Cicerone](https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/Cicerone_Data_Exploration/livetime_and_exposure.html).\n",
    "\n",
    "Here is the example of how to run [gtltcube](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtltcube.txt) (this step takes several minutes):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtltcube zmax=90\n",
    "    ./data/3C279_region_filtered_gti.fits\n",
    "    ./data/spacecraft.fits\n",
    "    ./data/3C279_ltcube.fits\n",
    "    0.025\n",
    "    1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">**NOTE**: Values such as 0.1 for `Step size in cos(theta)` are known to give unexpected reuslts. Use `0.09` instead."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " The livetime cube generated for this analysis can be found [here](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/Likelihood/3C279_ltcube.fits).\n",
    "\n",
    "Since [gtltcube](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtltcube.txt) produces a FITS file covering the entire sky, the output of this tool could be used for generating exposure maps for ROIs in other parts of the sky that have the same time interval selections.\n",
    "\n",
    "But use caution! The livetime cube MUST be regenerated if you change any part of the time interval selection. This can occur by changing the start or stop time of the events, or simply by changing the energy selection or zenith angle cut (as these produce a different set of GTIs from [gtmktime](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtmktime.txt)).\n",
    "\n",
    "See e.g. [Data Preparation](https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/Cicerone_Data_Exploration/Data_preparation.html) in the [Cicerone](https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/) or the [Data Preparation](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data_preparation.html) analysis thread.\n",
    "\n",
    "Although the [gtexpmap](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtexpmap.txt) application (see below) can generate exposure maps for Likelihood without a livetime file, using one affords a substantial time savings. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### b. Generate an exposure map.\n",
    "\n",
    "The tool [gtexpmap](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtexpmap.txt) creates an exposure map based on the event selection used on the input photon file and the livetime cube. The exposure map must be recalculated if the ROI, zenith, energy selection or the time interval selection of the events is changed. For more information about the exposure maps see the documentation in the Cicerone.\n",
    "\n",
    "Creating the exposure map using the [gtexpmap](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtexpmap.txt) tool, we have (this step can also take several minutes):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtexpmap\n",
    "    ./data/3C279_region_filtered_gti.fits\n",
    "    ./data/spacecraft.fits\n",
    "    ./data/3C279_ltcube.fits\n",
    "    ./data/3C279_unbin_expmap.fits\n",
    "    CALDB\n",
    "    30\n",
    "    120\n",
    "    120\n",
    "    20"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">**Note**: If trying to use `CALDB` in the example above does not work, just put the IRF name in explicitly. The appropriate IRF for this data set is `P8R3_SOURCE_V2`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we have chosen a 30 degree radius 'source region', while the acceptance cone radius specified for gtselect was 20 degrees. See the discussion of region selection in the [Cicerone](https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/Cicerone_Likelihood/Choose_Data.html). This is necessary to ensure that events from sources outside the ROI are accounted for owing to the size of the instrument point-spread function (PSF).\n",
    "\n",
    "Half-degree pixels are a nominal choice; smaller pixels should result in a more accurate evaluation of the diffuse source fluxes, but will also make the exposure map calculation itself lengthier.\n",
    "\n",
    "The number of energies specifies the number of logarithmically spaced intervals bounded by the energy range given in the DSS keywords. A general recommendation is 10 bins per decade. This is sufficient to accommodate the change in effective area with energy near 100 MeV.\n",
    "\n",
    "Here is one image plane of the exposure map we just created:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/images/Likelihood/3C279_expmap.png'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">**REMEMBER**: The exposure needs to be recalculated if the ROI, zenith angle, time, event class, or energy selections applied to the data are changed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Download the latest background models\n",
    "\n",
    "When you use the current Galactic diffuse emission model (`gll_iem_v07.fits`) in a likelihood analysis, you also want to use the corresponding model for the extragalactic isotropic diffuse emission, which includes the residual cosmic-ray background. The recommended isotropic model for point source analysis is `iso_P8R3_SOURCE_V2_v1.txt`.\n",
    "\n",
    "All the Pass 8 background models have been included in the Fermitools distribution in the `$(FERMI_DIR)/refdata/fermi/galdiffuse/` directory. If you use that path in your model, you should not have to download the diffuse models individually.\n",
    "\n",
    ">**NOTE**: Keep in mind that the isotropic model needs to agree with both the event class and event type selections you are using in your analysis. The iso_P8R3_SOURCE_V2_v1.txt isotropic spectrum is valid only for the latest response functions and only for data sets with front + back events combined. All of the most up-to-date background models along with a description of the models are available [here](https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The [gtlike](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtlike.txt) tool reads the source model from an XML file. The model file contains your best guess at the locations and spectral forms for the sources in your data. A source model can be created using the [model editor](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/modeleditor.txt) tool, by using the user contributed tool `make4FGLxml.py`, or by editing the file directly within a text editor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/analysis/user/make4FGLxml.py\n",
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/aux/4fgl/gll_iem_v07.fits\n",
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/aux/4fgl/iso_P8R3_SOURCE_V2_v1.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!mv make4FGLxml.py ./data/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given the dearth of bright sources in the extraction region we have selected, our source model file will be fairly simple, comprising only the Galactic and Extragalactic diffuse emission, and point sources to represent the blazars 3C 279 and 3C 273. Here is an example XML file:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```xml\n",
    "<?xml version=\"1.0\" ?>\n",
    "<source_library title=\"source library\">\n",
    "\n",
    "<!-- Point Sources -->\n",
    "\n",
    "<source name=\"3C 273\" type=\"PointSource\">\n",
    "<spectrum type=\"PowerLaw\">\n",
    "<parameter free=\"1\" max=\"1000.0\" min=\"0.001\" name=\"Prefactor\" scale=\"1e-09\" value=\"10\"/>\n",
    "<parameter free=\"1\" max=\"-1.0\" min=\"-5.0\" name=\"Index\" scale=\"1.0\" value=\"-2.1\"/>\n",
    "<parameter free=\"0\" max=\"2000.0\" min=\"30.0\" name=\"Scale\" scale=\"1.0\" value=\"100.0\"/>\n",
    "</spectrum>\n",
    "<spatialModel type=\"SkyDirFunction\">\n",
    "<parameter free=\"0\" max=\"360\" min=\"-360\" name=\"RA\" scale=\"1.0\" value=\"187.25\"/>\n",
    "<parameter free=\"0\" max=\"90\" min=\"-90\" name=\"DEC\" scale=\"1.0\" value=\"2.17\"/>\n",
    "</spatialModel>\n",
    "</source>\n",
    "\n",
    "<source name=\"3C 279\" type=\"PointSource\">\n",
    "<spectrum type=\"PowerLaw\">\n",
    "<parameter free=\"1\" max=\"1000.0\" min=\"0.001\" name=\"Prefactor\" scale=\"1e-09\" value=\"10\"/>\n",
    "<parameter free=\"1\" max=\"-1.0\" min=\"-5.0\" name=\"Index\" scale=\"1.0\" value=\"-2\"/>\n",
    "<parameter free=\"0\" max=\"2000.0\" min=\"30.0\" name=\"Scale\" scale=\"1.0\" value=\"100.0\"/>\n",
    "</spectrum>\n",
    "<spatialModel type=\"SkyDirFunction\">\n",
    "<parameter free=\"0\" max=\"360\" min=\"-360\" name=\"RA\" scale=\"1.0\" value=\"193.98\"/>\n",
    "<parameter free=\"0\" max=\"90\" min=\"-90\" name=\"DEC\" scale=\"1.0\" value=\"-5.82\"/>\n",
    "</spatialModel>\n",
    "</source>\n",
    "\n",
    "<!-- Diffuse Sources -->\n",
    "<source name=\"gll_iem_v07\" type=\"DiffuseSource\">\n",
    "<spectrum type=\"PowerLaw\">\n",
    "<parameter free=\"1\" max=\"10\" min=\"0\" name=\"Prefactor\" scale=\"1\" value=\"1\"/>\n",
    "<parameter free=\"0\" max=\"1\" min=\"-1\" name=\"Index\" scale=\"1.0\" value=\"0\"/>\n",
    "<parameter free=\"0\" max=\"2e2\" min=\"5e1\" name=\"Scale\" scale=\"1.0\" value=\"1e2\"/>\n",
    "</spectrum>\n",
    "<spatialModel file=\"./gll_iem_v07.fits\" type=\"MapCubeFunction\">\n",
    "<parameter free=\"0\" max=\"1e3\" min=\"1e-3\" name=\"Normalization\" scale=\"1.0\" value=\"1.0\"/>\n",
    "</spatialModel>\n",
    "</source>\n",
    "<source name=\"iso_P8R3_SOURCE_V2_v1\" type=\"DiffuseSource\">\n",
    "<spectrum file=\"./iso_P8R3_SOURCE_V2_v1.txt\" type=\"FileFunction\" apply_edisp=\"false\">\n",
    "<parameter free=\"1\" max=\"10\" min=\"1e-2\" name=\"Normalization\" scale=\"1\" value=\"1\"/>\n",
    "</spectrum>\n",
    "<spatialModel type=\"ConstantValue\">\n",
    "<parameter free=\"0\" max=\"10.0\" min=\"0.0\" name=\"Value\" scale=\"1.0\" value=\"1.0\"/>\n",
    "</spatialModel>\n",
    "</source>\n",
    "</source_library>\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The XML file used for this example is `3C279input_model.xml` (in the cell below). For more details on the available XML models, see:\n",
    "* Descriptions of available [Spectral and Spatial Models](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/source_models.html)\n",
    "* Examples of [XML Model Definitions for Likelihood](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/xml_model_defs.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/Likelihood/3C279input_model.xml"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!mv 3C279input_model.xml ./data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### XML for Extended Sources\n",
    "\n",
    "If your region includes any extended sources, you will need to add one or more extended sources to your XML model.\n",
    "\n",
    "Download the extended source templates from the [LAT Catalog](https://fermi.gsfc.nasa.gov/ssc/data/access/lat/4yr_catalog/) page (look for \"Extended Source template archive\"). Extract the archive in the directory of your choice and note the path to the template files, which have names like `W44.fits` and `VelaX.fits`. You will need to modify your XML model so that the path to the template file is correct.\n",
    "\n",
    "In addition, you will need to add the modifier `map_based_integral='true'` to the `<spatialModel>` atttribute. This modifier tells [gtlike](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtlike.txt) to use the proper integration method when evaluating the extended source.\n",
    "\n",
    "Here is an example of the proper format for an extended source XML entry for Unbinned Likelihood analysis:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```xml\n",
    "<?xml version=\"1.0\" standalone=\"no\"?> \n",
    "<source_library title=\"source library\"> \n",
    "\n",
    "<source name=\"W44\" type=\"DiffuseSource\"> \n",
    "<spectrum type=\"LogParabola\" normPar=\"norm\"> \n",
    "<parameter free=\"1\" max=\"100000\" min=\"1e-05\" name=\"norm\" scale=\"1e-11\" value=\"2.891247798\"/> \n",
    "<parameter free=\"1\" max=\"5\" min=\"0\" name=\"alpha\" scale=\"1\" value=\"2.399924084\"/> \n",
    "<parameter free=\"1\" max=\"5\" min=\"-1\" name=\"beta\" scale=\"1\" value=\"0.2532196784\"/> \n",
    "<parameter free=\"0\" max=\"300000\" min=\"20\" name=\"Eb\" scale=\"1\" value=\"1737.352167\"/> \n",
    "</spectrum> \n",
    "<spatialModel file=\"./W44.fits\" map_based_integral=\"true\" type=\"SpatialMap\"> \n",
    "<parameter free=\"0\" max=\"1000\" min=\"0.001\" name=\"Prefactor\" scale=\"1\" value=\"1\"/> \n",
    "</spatialModel> \n",
    "</source> \n",
    "\n",
    "</source_library> \n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 6. Compute the diffuse source responses\n",
    "\n",
    "The diffuse source responses are computed by the [gtdiffrsp](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtdiffrsp.txt) tool.\n",
    "\n",
    "The source model XML file must contain all of the diffuse sources to be fit. [gtdiffrsp](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtdiffrsp.txt) will add one column to the event data file for each diffuse source. The diffuse response depends on the instrument response function (IRF), which must be in agreement with the selection of events, i.e. the event class and event type we are using in our analysis. In this case we are using the standard SOURCE class (`evclass=128`) and event type 3 (`FRONT+BACK`), so we should use the IRF `P8R3_SOURCE_V2` (which is the most commonly used in standard analysis).\n",
    "\n",
    "Alternatively, you can just use `CALDB` for the IRF and the appropriate IRF will be automatically applied according to the previous event selection you performed during the previous steps of the analysis.\n",
    "\n",
    "If the diffuse responses are not precomputed using [gtdiffrsp](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtdiffrsp.txt), then the [gtlike](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtlike.txt) tool will compute them at runtime (during the next step). However, as this step is very computationally intensive (often taking ~hours to complete), and it is very likely you will need to run [gtlike](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtlike.txt) more than once, it is probably wise to precompute these quantities.\n",
    "\n",
    "Here is an example of running [gtdiffrsp](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtdiffrsp.txt):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtdiffrsp\n",
    "    ./data/3C279_region_filtered_gti.fits\n",
    "    ./data/spacecraft.fits\n",
    "    ./data/3C279input_model.xml\n",
    "    CALDB"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">**Note**: If trying to use CALDB in the example above does not work, just put the IRF name in explicitly. The appropriate IRF for this data set is `P8R3_SOURCE_V2`.\n",
    "\n",
    "There are a number of IRFs distributed with the Fermitools. Some of these IRFs, like `P8R3_SOURCE_V2::FRONT`, are designed to address only a subset of the events in the dataset, and thus would require you to make additional cuts before they could be used.\n",
    "\n",
    "Once you have run [gtdiffrsp](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtdiffrsp.txt), you can find the current names of the diffuse response column using the **FTOOLS** with:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```bash\n",
    "fkeyprint <event file name>.fits DIFRSP\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "fkeyprint\n",
    "    ./data/3C279_region_filtered_gti.fits\n",
    "    DIFRSP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first part tells you the IRF name and, after the double underscore, the name of the diffuse component associated with the response.\n",
    "\n",
    "If you are using a different IRF or diffuse model, you MUST compute a custom diffuse response prior to continuing the analysis.\n",
    "\n",
    "Now we are ready to perform a likelihood analysis. You can get the resulting file, including the diffuse responses, [here](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/Likelihood/3C279_events_gti_diffs.fits)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 7. Run gtlike\n",
    "\n",
    "We are now ready to run the [gtlike](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtlike.txt) application.\n",
    "\n",
    "You may need the results of the likelihood fit to be output in XML model form (e.g. to use in generating a test statistic map).\n",
    "\n",
    "To obtain an XML model output in addition to the standard results files, use the `sfile` parameter on the command line (as shown below) to designate the output XML model filename."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtlike refit=yes plot=yes sfile=3C279output_model.xml\n",
    "    UNBINNED\n",
    "    ./data/spacecraft.fits\n",
    "    ./data/3C279_region_filtered_gti.fits\n",
    "    ./data/3C279_unbin_expmap.fits\n",
    "    ./data/3C279_ltcube.fits\n",
    "    ./data/3C279input_model.xml\n",
    "    CALDB\n",
    "    NEWMINUIT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most of the entries prompted for are fairly obvious.\n",
    "\n",
    "In addition to the various XML and FITS files, the user is prompted for a choice of IRFs, the type of statistic to use, the optimizer, and some output file names.\n",
    "\n",
    "The statistics available are:\n",
    "* **UNBINNED**: This is a standard unbinned analysis, described in this tutorial, to be used for short timescale or low source count data. If this option is chosen then parameters for the spacecraft file, event file, and exposure file must be given.\n",
    "\n",
    "* **BINNED**: This analysis is used for long timescale or high-density data (such as in the Galactic plane) which can cause memory errors in the unbinned analysis. See an explanation in the [Binned Likelihood Tutorial](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/binned_likelihood_tutorial.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are five optimizers from which to choose: `DRMNGB`, `DRMNFB`, `NEWMINUIT`, `MINUIT` and `LBFGS`.\n",
    "\n",
    "Generally speaking, the faster way to find the parameter estimates is to use `DRMNGB` (or `DRMNFB`) to find initial values and then use `MINUIT` (or `NEWMINUIT`) to find more accurate results.\n",
    "\n",
    "If you have trouble achieving convergence at first, you can loosen your tolerance by setting the hidden parameter `ftol` on the command line. (The default value for `ftol` is 0.01.)\n",
    "\n",
    "The application proceeds by reading in the spacecraft and event data, and if necessary, computing event responses for each diffuse source.\n",
    "If you get an error similar to this:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "Caught St11range_error at the top level: Requested energy, 911163, \n",
    "lies outside the range of the input file, 34.171, 877938\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This means that we are trying to do our analysis beyond the maximum energy covered by the data sample. We should have checked that it was not the case during the data selection (as recommended in step 1).\n",
    "\n",
    "To solve this issue we need to apply a harder cut in our event selection, setting the maximum energy within the energy range of available photons in your data sample.\n",
    "\n",
    "Fortunately, this new event selection will not affect the time calculation (previous steps do not need to be repeated).\n",
    "\n",
    "Here is the output from our fit with [gtlike](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtlike.txt):"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "Minuit did successfully converge. \n",
    "\n",
    "(MUCH OUTPUT OMITTED.) \n",
    "\n",
    "Computing TS values for each source (4 total) \n",
    "....! \n",
    "\n",
    "Photon fluxes are computed for the energy range 100 to 500000 MeV \n",
    "\n",
    "3C 273: \n",
    "Prefactor: 9.53263 +/- 0.262701 \n",
    "Index: -2.63212 +/- 0.0212044 \n",
    "Scale: 100 \n",
    "Npred: 6369.65 \n",
    "ROI distance: 10.4409 \n",
    "TS value: 7810.96 \n",
    "Flux: 5.85022e-07 +/- 1.14156e-08 photons/cm^2/s \n",
    "\n",
    "3C 279: \n",
    "Prefactor: 7.3524 +/- 0.189805 \n",
    "Index: -2.20869 +/- 0.0144918 \n",
    "Scale: 100 \n",
    "Npred: 7420.98 \n",
    "ROI distance: 0 \n",
    "TS value: 13068.8 \n",
    "Flux: 6.08821e-07 +/- 1.08338e-08 photons/cm^2/s \n",
    "\n",
    "gll_iem_v07: \n",
    "Prefactor: 1.20506 +/- 0.0160704 \n",
    "Index: 0 \n",
    "Scale: 100 \n",
    "Npred: 95322.2 \n",
    "Flux: 0.000627754 +/- 8.37036e-06 photons/cm^2/s \n",
    "\n",
    "iso_P8R3_SOURCE_V2: \n",
    "Normalization: 1.05663 +/- 0.0281431 \n",
    "Npred: 47447.5 \n",
    "Flux: 0.000129821 +/- 3.45716e-06 photons/cm^2/s \n",
    "\n",
    "WARNING: Fit may be bad in range [100, 549.28] (MeV) \n",
    "WARNING: Fit may be bad in range [213340, 326604] (MeV) \n",
    "\n",
    "Total number of observed counts: 156560 \n",
    "Total number of model events: 156560 \n",
    "\n",
    "-log(Likelihood): 1670346.65 \n",
    "\n",
    "Writing fitted model to 3C279output_model.xml\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we selected `plot=yes` in the command line, a plot of the fitted data appears.\n",
    "\n",
    "<img src='https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/images/Likelihood/3C279_spectral_fit.png'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the first plot, the counts/MeV vs MeV are plotted. The points are the data, and the lines are the models. Error bars on the points represent sqrt(Nobs) in that band, where Nobs is the observed number of counts. The black line is the sum of the models for all sources. \n",
    "\n",
    "The colored lines follow the sources as follows:\n",
    "\n",
    "* Black - summed model\n",
    "* Red - first source in results.dat file (see below\n",
    "* Green - second source\n",
    "* Blue - third source\n",
    "* Magenta - fourth source\n",
    "\n",
    "In our case the 3C273 is the one in red, 3C279 is the one in green, the extragalactic background is blue, and the galactic background is magenta. If you have more sources, the colors are reused in the same order."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second plot gives the residuals between your model and the data. Error bars here represent (sqrt(Nopbs))/Npred, where Npred is the predicted number of counts in each band based on the fitted model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To assess the quality of the fit, look first for the words at the top of the output `<Optimizer> did successfully converge.` Successful convergence is a minimum requirement for a good fit.\n",
    "\n",
    "Next, look at the energy ranges that are generating warnings of bad fits. If any of these ranges affect your source of interest, you may need to revise the source model and refit. You can also look at the residuals on the plot (bottom panel). If the residuals indicate a poor fit overall, you should consider changing your model file, perhaps by using a different source model definition and/or adding new sources in the ROI, and then refit the data.\n",
    "\n",
    "If the fits and spectral shapes are good, but could be improved, you may wish to simply update your model file to hold some of the spectral parameters fixed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, by fixing the spectral model for 3C 273, you may get a better quality fit for 3C 279. Close the plot and you will be asked if you wish to refit the data:\n",
    "```\n",
    "Refit? [y]n \n",
    "Elapsed CPU time: 111.819315\n",
    "prompt>\n",
    "```\n",
    "Here, hitting \"return\" will instruct the application to fit again. We are happy with the result, so we type `n` and end the fit."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Results\n",
    "\n",
    "When it completes, [gtlike](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtlike.txt) generates a standard output XML file with the results of the fit.\n",
    "\n",
    "If you re-run the tool in the same directory, these files will be overwritten by default. Use the `clobber=no` option on the command line to keep from overwriting the output files.\n",
    "\n",
    "Unfortunately, the fit details and the value for the -log(likelihood) are not recorded in the automatic output files. You should consider logging the output to a text file for your records by using `gtlike ... > fit_data.txt` (or something similar) with your [gtlike](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtlike.txt) command.\n",
    "\n",
    "Be aware, however, that this will make it impossible to request a refit when the likelihood process completes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Example:\n",
    "```bash\n",
    "gtlike plot=yes sfile=3C279output_model.xml > fit_data.txt\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, we used the `sfile` parameter to request that the model results be written to an output XML file, `3C279output_model.xml`.\n",
    "\n",
    "This file contains the source model results that were written to `results.dat` at the completion of the fit.\n",
    "\n",
    ">**NOTE**: If you have specified an output XML model file and you wish to modify your model while waiting at the `Refit? [y]` prompt, you will need to copy the results of the output model file to your input model before making those modifications.\n",
    "\n",
    "Errors reported with each value in the results.dat file are 1Ïƒ estimates (based on the inverse-Hessian at the optimum of the log-likelihood surface)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Other Useful Hidden Parameters\n",
    "\n",
    "If you are scripting and wish to generate multiple output files without overwriting, the `results` and `specfile` parameters allow you to specify output filenames for the `results.dat` and `counts_spectra.fits` files respectively.\n",
    "\n",
    "If you do not specify a source model output file with the `sfile` parameter, then the input model file will be overwritten with the latest fit. This is convenient as it allows the user to edit that file while the application is waiting at the `Refit? [y]` prompt so that parameters can be adjusted and set free or fixed. This would be similar to the use of the `newpar`, `freeze`, and `thaw` commands in [XSPEC](http://heasarc.gsfc.nasa.gov/docs/xanadu/xspec/index.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 8. Make test-statistic maps\n",
    "\n",
    "Ultimately, one would like to find sources near the detection limit of the instrument.\n",
    "\n",
    "To do this, you model the strongest, most obvious sources (with some theoretical prejudice as to the true source positions, e.g., assuming that most variable high Galactic latitude sources are blazars which can be localized by radio, optical, or X-ray observations), and then create `Test-statistic maps` (TS maps) to search for unmodeled point sources.\n",
    "\n",
    "These TS maps are created by moving a putative point source through a grid of locations on the sky and maximizing -log(likelihood) at each grid point, with the other, stronger, and presumably well-identified sources included in each fit. New, fainter sources are then identified at local maxima of the TS map."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The [gttsmap](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gttsmap.txt) tool can be run with or without an input source model. However, for a useful visualization of the results of the fit, it is recommended you use the output model file from [gtlike](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtlike.txt). The file must be edited so all parameters are fixed (by setting the `free` attribute to 0 for each parameter) otherwise [gttsmap](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gttsmap.txt) will attempt a refit of the entire model at every point on the grid.\n",
    "\n",
    "To see the field with the fitted sources _removed_ (i.e. a residuals map), fix all point source parameters before running the TS map. See residuals map wget cell below.\n",
    "\n",
    "To see the field with the fitted sources _included_, edit the model to remove all but the diffuse components. See sources map wget cell below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# For residuals map: 3C279output_model_resid.xml\n",
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/Likelihood/3C279output_model_resid.xml"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# For sources TS map: 3C279output_model_diffonly.xml\n",
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/Likelihood/3C279output_model_diffonly.xml"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!mv *.xml ./data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In both cases, leave the Galactic diffuse prefactor (called `Value` in the model file) and the isotropic diffuse normalization parameters free during the fit.\n",
    "\n",
    "Running [gttsmap](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gttsmap.txt) is **extremely time-consuming**, as the tool is performing a likelihood fit for all events at every pixel position.\n",
    "\n",
    "One way to reduce the time required for this step is to use very coarse binning and/or a very small region.\n",
    "\n",
    "In the following example, we run a TS map for the central 20x20 degree region of our data file with .25 degree bins. This results in 6400 maximum likelihood calculations.\n",
    "\n",
    "The run time for each of the maps discussed below was approximately one week.\n",
    "\n",
    "Here is an example of how to run the [gttsmap](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gttsmap.txt) tool to look for additional sources:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gttsmap\n",
    "    UNBINNED\n",
    "    ./data/spacecraft.fits\n",
    "    ./data/3C279_region_filtered_gti.fits\n",
    "    ./data/3C279_unbin_expmap.fits\n",
    "    ./data/3C279_ltcube.fits\n",
    "    ./data/3C279output_model_resid.xml\n",
    "    CALDB\n",
    "    NEWMINUIT\n",
    "    ./data/3C279_tsmap_resid.fits\n",
    "    80\n",
    "    80\n",
    "    0.25\n",
    "    CEL\n",
    "    193.98\n",
    "    -5.82\n",
    "    AIT\n",
    "\n",
    "#### Parameters\n",
    "#### Takes a long time to run!\n",
    "# Statistic to use (BINNED|UNBINNED)\n",
    "# Spacecraft file\n",
    "# Event file \n",
    "# Unbinned exposure map\n",
    "# Exposure hypercube file\n",
    "# Source model file\n",
    "# Response functions to use\n",
    "# Optimizer (DRMNFB|NEWMINUIT|MINUIT|DRMNGB|LBFGS)\n",
    "# TS map file name\n",
    "# Number of X axis pixels\n",
    "# Number of Y axis pixels\n",
    "# Image scale (in degrees/pixel)\n",
    "# Coordinate system (CEL|GAL)\n",
    "# X-coordinate of image center in degrees (RA or l)\n",
    "# Y-coordinate of image center in degrees (Dec or b)\n",
    "# Projection method (AIT|ARC|CAR|GLS|MER|NCP|SIN|STG|TAN)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">**Note**: If trying to use CALDB in the example above does not work, just put the IRF name in explicitly. The appropriate IRF for this data set is `P8R3_SOURCE_V2`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because generating TS maps takes a long time, you may wish to download both the residuals and source files in the cells below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Residual file\n",
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/Likelihood/3C279_tsmap_resid.fits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Source file\n",
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/Likelihood/3C279_tsmap_source.fits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!mv *.fits ./data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output from the fit is below.\n",
    "\n",
    "In the left panel, only diffuse sources were included in the analysis.\n",
    "\n",
    "The right panel shows the same field, but with the point sources (3C 279 and 3C 273) included in the model, and thus not included in the output image.\n",
    "\n",
    "This gives a pseudo-residuals map.\n",
    "\n",
    "The bright source clearly seen in the lower part of this image is a recently identified millisecond pulsar. The location of the cursor in the image indicates the TS value at each position.\n",
    "\n",
    "<img src='https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/images/Likelihood/3C279_tsmap_tiled.png'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">**NOTE**: There may be a number of black pixels in these images where the likelihood fit did not converge. You should try adjusting the tolerance or using a different minimizer to keep this from happening."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, the data set was small enough that an **unbinned likelihood** analysis is possible.\n",
    "\n",
    "With longer data sets, however, you may run into memory allocation errors, at which time it is necessary to move to **binned analysis** for your region. A typical memory allocation error looks like:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "gtlike(10869) malloc: *** mmap(size=2097152) failed (error code=12)*** error: can't allocate region*** set a breakpoint in malloc_error_break to debug\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you encounter this error, then [Binned Likelihood](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/binned_likelihood_tutorial.html) analysis is for you!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
